function [hfile, densm] = srf_Morph(hfile, niter, force, type, opts)
% SRF::Morph  - apply vertex morphing
%
% FORMAT:       [srf, densm] = srf.Morph([niter, force [, type [, opts]])
%
% Input fields:
%
%       niter       number of iterations (default 1)
%       force       smoothing force applied to morphing (default: 0.07)
%       type        either of 'even', {'smooth'}
%       opts        optional settings
%        .areac     if given and true, keep area constant
%        .distc     additional distortion correction force (>0 ... 3)
%        .pbar      progress bar flag (default: from ini)
%        .sphere    additionally applied to-sphere force
%        .title     display title for progress bar
%
% Output fields:
%
%       srf         altered object
%       densm       density map SMP
%
% Note: this method requires the MEX file meshmorph to be compiled!

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% global configuration
global bvqxfile_config;
stc = bvqxfile_config.settings;
if ~isfield(stc, 'Morphing')
    stc = struct('Morphing', struct);
end
stc = stc.Morphing;
if isfield(stc, 'ShowProgress')
    spbar = stc.ShowProgress;
else
    spbar = true;
end
if isfield(stc, 'ShowSRF')
    spsrf = stc.ShowSRF;
else
    spsrf = false;
end

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
if nargin < 5 || ...
    numel(opts) ~= 1 || ...
   ~isstruct(opts)
    opts = struct;
end
if nargin < 2 || ...
    numel(niter) ~= 1 || ...
   ~isa(niter, 'double') || ...
	isnan(niter) || ...
    niter < 1 || ...
    niter > 1e5
    niter = 1;
else
    niter = floor(niter);
end
if nargin < 3 || ...
    numel(force) ~= 1 || ...
   ~isa(force, 'double') || ...
	isnan(force) || ...
    force < -16 || ...
    force > 16 || ...
    force == 0
    force = 0.07;
end
opts.force = force;
if nargin < 4 || ...
    isempty(type) || ...
   ~ischar(type) || ...
   ~any(strcmpi(type(:)', {'even', 'smooth'}))
    opts.type = 'smooth';
else
    opts.type = lower(type(:)');
end
if isfield(opts, 'areac')
    if numel(opts.areac) == 1
        if (islogical(opts.areac) && ...
            opts.areac) || ...
           (isnumeric(opts.areac) && ...
            opts.areac ~= 0)
            opts.areac = double(0 + (opts.areac && true));
        end
    else
        opts.areac = 1;
    end
end
if isfield(opts, 'pbar') && ...
    islogical(opts.pbar) && ...
    numel(opts.pbar) == 1
    spbar = opts.pbar;
end
if ~isfield(opts, 'title') || ...
   ~ischar(opts.title) || ...
    isempty(opts.title)
    opts.title = [upper(opts.type(1)) opts.type(2:end) 'ing'];
end

% get content
sc = bvqxfile_getscont(hfile.L);
bc = sc.C;
c = bc.VertexCoordinate;
n = bc.Neighbors;
[fn{1:3}] = fileparts(sc.F);
fn{2} = [fn{2} fn{3}];
if numel(fn{2}) > 28
    fn{2} = [fn{2}(1:13) '...' fn{2}(end-12:end)];
end

% determine progress bar capabilities
pbar = [];
if spbar
    try
        pbar = BVQXprogress;
        BVQXprogress(pbar, 'setposition', [80, 120, 640, 36]);
        BVQXprogress(pbar, 'settitle', sprintf('%s - morphing', fn{2}));
        BVQXprogress(pbar, 0, ...
            sprintf('%s (%d steps)...', opts.title, niter), 'visible', 0, niter + 40);
    catch
        % do nothing
    end
end

% also show SRF
if spsrf
    srf_Show(hfile);
    bc = bvqxfile_getcont(hfile.L);
end

% if pbar, loop while niter >= stepsize
stepsize = min(250, max(50, ceil(niter ^ (2 / 3))));
done = 0;
if ~isempty(pbar) || ...
    spsrf
    while niter >= stepsize
        opts.niter = stepsize;
        try
            c = mesh_morph(c, n, bc.TriangleVertex, opts);
        catch
            closebar(pbar);
            error( ...
                'BVQXfile:MEXError', ...
                'Error morphing mesh (after %d pre-iterations): %s.', ...
                done, lasterr ...
            );
        end
        niter = niter - stepsize;
        done = done + stepsize;
        if ~isempty(pbar)
            BVQXprogress(pbar, done, ...
              sprintf('%s (%d of %d steps)...', opts.title, done, done + niter));
        end
        if spsrf
            bc.VertexCoordinate = c;
            bvqxfile_setcont(hfile.L, bc);
            srf_RecalcNormals(hfile);
            srf_Show(hfile);
        end
    end
end
opts.niter = niter;

% last call to morphmesh
try
    if nargout > 1
        [c, densmv] = mesh_morph(c, n, bc.TriangleVertex, opts);
    else
        c = mesh_morph(c, n, bc.TriangleVertex, opts);
    end
catch
    if ~isempty(pbar)
        closebar(pbar);
    end
    error( ...
        'BVQXfile:MEXError', ...
        'Error morphing mesh (after %d pre-iterations): %s.', ...
        done, lasterr ...
    );
end

% put back
bc.VertexCoordinate = c;
bvqxfile_setcont(hfile.L, bc);

% progress bar
if ~isempty(pbar)
    BVQXprogress(pbar, niter + done, ...
        sprintf('%s (re-calculcating normals)', opts.title));
end

% and recalc normals after morphing
srf_RecalcNormals(hfile);

% delete figure with surface
if spsrf
    h = aft_Handles(hfile);
    if isfield(h, 'Surface') && ...
        numel(h.Surface) == 1 && ...
        ishandle(h.Surface)
        h = h.Surface;
        while ~strcmpi(get(h, 'Type'), 'figure')
            h = get(h, 'Parent');
            if strcmpi(get(h, 'Type'), 'root')
                warning( ...
                    'BVQXtools:InternalError', ...
                    'Invalid UI handle structure. Leaving surface.' ...
                );
                break;
            end
        end
        if strcmpi(get(h, 'Type'), 'figure')
            delete(h);
            aft_DeleteHandle(hfile, 'Surface');
        end
    end
end

% produce true output map
if nargout > 1
    
    % progress bar
    if ~isempty(pbar)
        BVQXprogress(pbar, niter + done + 20, ...
            sprintf('%s (creating density SMP)', opts.title));
    end

    % create map
    densmv = densmv.DensitySMPData;
    md = mean(densmv);
    sd = std(densmv);
    densm = BVQXfile('new:smp');
    dbc = bvqxfile_getcont(densm.L);
    dbc.NrOfVertices = numel(densmv);
    dbc.NameOfOriginalSRF = sc.F;
    dbc.Map.LowerThreshold = max(0.01, md - 1.5 * sd);
    dbc.Map.UpperThreshold = md + 1.5 * sd;
    dbc.Map.UseValuesAboveThresh = 0;
    dbc.Map.DF1 = 6;
    dbc.Map.BonferroniValue = numel(densmv);
    dbc.Map.Name = 'Mesh density map (average area per vertex)';
    dbc.Map.SMPData = densmv;
    dv = (densmv - md) / sd;
    dbc.Map(2) = dbc.Map(1);
    dbc.Map(2).LowerThreshold = 1;
    dbc.Map(2).UpperThreshold = 2.5;
    dbc.Map(2).UseValuesAboveThresh = 1;
    dbc.Map(2).DF1 = 6;
    dbc.Map(2).BonferroniValue = numel(densmv);
    dbc.Map(2).Name = 'Mesh density map (normalized average area per vertex)';
    dbc.Map(2).SMPData = dv;
    bvqxfile_setcont(densm.L, dbc);
end

% clean up progress bar
if ~isempty(pbar)
    closebar(pbar);
end
